% The following are a set of test problems for the simple discrete
% ordinates program "sn".  They consist of both fixed source and eigenvalue
% problem.  They are in brief:
%
% Problems 1, 2, and 3
%   Slab reactor benchmarks from Ilas and Rahnema, "A Heterogeneous Coarse
%   Mesh Transport Method", Trans. Th. Stat. Phys., 32, 445 - 471, (2003) 
%   They also appear in S. Mosher's PhD thesis.  A total of four assembly
%   types, each containing six material regions are used to construct three
%   different core configurations with increasing heterogeneity.  The
%   problems are two-group.
%
%   The benchmark keff's are taken from Mosher.  A uniform fine mesh of 
%   0.579 cm was used with an S32 Gauss-Legendre quadrature.  Fluxes are
%   converged to less than 1e-6.  Inners were limited to 20 per group per
%   outer iteration.  
%                          k      total inners    total outers
%       Problem 1       1.258247      1825            57
%       Problem 2       1.007066      1943            58
%       Problem 3       0.805372      2424            74
%   Our test code gets these values within convergence for the same
%   discretization.
%
%   Note, for eigenproblems, the "src" variable is used to make a fission
%   density guess at the coarse mesh level.  Just use the first group (as
%   it's normalized).
%
% Problems 4, 5, 6, 7, and 8
%   A one-group, homogeneous, highly-scattering slab.  For various coarse
%   mesh size and fixed fine mesh size, the efficacy (and stability) of
%   CMR, CMFD, and CMSC can be identified.  Used S8 and eps = 1e-6 for
%   source.
%                                 
%   Problem |    SI  |  CMR  | CMFD  | CMSC  | mesh(mfp)
%      4    |   910  |   67  |diverge|  74   | 5
%      5    |   910  |   29  |diverge|  37   | 2.5   
%      6    |   910  |   19  |    17 |  21   | 1.25   (stable up to ~1.67mfp)
%      7    |   910  |  186  |    10 |  38   | 0.625   
%      8    |   910  |diverge|     9 |diverge| 0.3125 
%
% Strangely, CMSC diverges for Problem 8 using S8 but not for S2, S16, or
% S32.
%
% Note, for 0.625 mfp mesh and c=0.9999, CMFD still takes just 7 its for the
% 1e-6 convergence on phi whereas unaccelerated run takes 14441 its.
%

p = 12; % PROBLEM IDENTIFIER

if ( p<4 || p==9 )
    base   = [ 1.1580 4.4790 7.8000 11.1210 14.4420 15.6000 ]; % an "assembly"
    xcm    = [ 0.0  base  base+15.6  base+15.6*2 base+15.6*3 base+15.6*4 base+15.6*5 base+15.6*6 ];
    basef  = [ 1 2 2 2 2 1 ]*2; % fine mesh discretization of assembly
    xfm    = [ basef basef basef basef basef basef basef ];
    bases  = [ 0 1 1 1 1 0
               0 0 0 0 0 0 ];
    src    = [ bases bases bases bases bases bases bases ];
    mtA    = [ 1 2 3 3 2 1 ];
    mtB    = [ 1 2 2 2 2 1 ];
    mtC    = [ 1 2 4 4 2 1 ];
    mtD    = [ 1 4 4 4 4 1 ];
    data   = [  0.1890 0.0000 0.0000 0.0000 0.1507 0.0380 % water
                1.4633 0.0000 0.0000 0.0000 0.0000 1.4536
                0.2263 0.0000 0.0067 1.0000 0.2006 0.0161 % fuel I
                1.0119 0.0000 0.1241 0.0000 0.0000 0.9355
                0.2252 0.0000 0.0078 1.0000 0.1995 0.0156 % fuel II
                0.9915 0.0000 0.1542 0.0000 0.0000 0.9014
                0.2173 0.0000 0.0056 1.0000 0.1902 0.0136 % fuel III
                1.0606 0.0000 0.0187 0.0000 0.0000 0.5733
             ] ;
    numm = 4; numg = 2; ord = 32; ptype = 1; 
elseif ( p < 9 )
    data   = [ 1.00   0.50  0 0  0.50    % mid scatter
               1.50   1.50  0 0  0.00    % pure absorber
               0.50   0.30  0 0  0.20    
               1.00   0.01  0 0  0.99    % high scatter
               1.00   0.00  0 0  1.0  ]; % pure scatter  
    numm = 5; numg = 1; ord = 2; ptype = 0; 
end

% EIGENVALUE EXAMPLES
if ( p==1 )
    % configuration 1:  A B A B A B A,
    mt     = [ mtA mtB mtA mtB mtA mtB mtA ];
elseif ( p==2 )
    % configuration 2:  A C A C A C A
    mt     = [ mtA mtC mtA mtC mtA mtC mtA ];
elseif ( p==3) 
    % configuration 3:  A B A B A B A
    mt     = [ mtA mtD mtA mtD mtA mtD mtA ];
% FIXED SOURCE EXAMPLES    
elseif (p>3 && p<9) 
    if ( p == 4 )
        xcm    = 0:5:100;
        xfm    = 32*ones(1,length(xcm)-1);
    elseif ( p == 5 )
        xcm    = 0:2.5:100;
        xfm    = 16*ones(1,length(xcm)-1);
    elseif ( p == 6 )
        xcm    = 0:1.25:100;
        xfm    = 8*ones(1,length(xcm)-1);
    elseif ( p == 7 )
        xcm    = 0:0.625:100;
        xfm    = 4*ones(1,length(xcm)-1);
    else
        xcm    = 0:0.3125:100;
        xfm    = 2*ones(1,length(xcm)-1);
    end
    src    = ones(1,length(xcm)-1);
    mt     = 4*ones(1,length(xcm)-1);
elseif ( p == 9 )
    % configuration 1:  A B A B A B A (neglecting fission)
    mt     = [ mtA mtB mtA mtB mtA mtB mtA ];
    numm = 4; numg = 2; ord = 32; ptype = 0; 
elseif ( p == 10 )
    clc
    % simple 1-d eigenvalue.  keff=1.1873965(9),138/8 outer w/,w/o cmd
    xcm  = 0:0.0625:20;
    xfm  = 4*ones(1,length(xcm)-1);
    src  = ones(1,length(xcm)-1);
    mt   = ones(1,length(xcm)-1);
    data = [ 1.40 0.50 0.60 1.0  0.90 ]; % kin = nuSigF/SigA = 1.2
    ptype = 1; ord = 32; numg =1; numm=1;
elseif ( p == 11 )
    clc
    % simple 1-d eigenvalue.  keff=1.1873965(9),138/8 outer w/,w/o cmd
    d = 0.03125;
    xcm  = 0:d:40;
    xfm  = 2*ones(1,length(xcm)-1);
    src  = ones(1,length(xcm)-1);
    tmp  = (length(xcm)-1)/16;
    
    mt   = [ 4*ones(1,tmp) 1*ones(1,tmp) 4*ones(1,tmp) 2*ones(1,tmp) ...
             4*ones(1,tmp) 3*ones(1,tmp) 4*ones(1,tmp) 1*ones(1,tmp) ...
             4*ones(1,tmp) 1*ones(1,tmp) 4*ones(1,tmp) 2*ones(1,tmp) ...
             4*ones(1,tmp) 3*ones(1,tmp) 4*ones(1,tmp) 1*ones(1,tmp)];
%     src   =[ 0*ones(1,tmp) 0*ones(1,tmp) 0*ones(1,tmp) 1*ones(1,tmp) ...
%              0*ones(1,tmp) 0*ones(1,tmp) 0*ones(1,tmp) 0*ones(1,tmp) ...
%              0*ones(1,tmp) 0*ones(1,tmp) 0*ones(1,tmp) 1*ones(1,tmp) ...
%              0*ones(1,tmp) 0*ones(1,tmp) 0*ones(1,tmp) 0*ones(1,tmp)];         
    data = [ 1.00 0.50 0.60 1.0  0.50
             1.00 0.50 0.55 1.0  0.50 
             1.00 0.50 0.55 1.0  0.50 
             1.00 0.50 0.63 1.0  0.50     ]; % kin = nuSigF/SigA = 1.2    
    ptype = 1; ord = 2; numg =1; numm=4;    
    % FINAL keff = 1.20999188 in 452 iterations
    % FINAL keff = 1.20999190 in 5 iterations
elseif ( p == 12 )
    xcm = 0:0.5:20;
    xfm = 100*ones(1,length(xcm)-1); nc2 = (length(xcm)-1)/2;
    src = zeros(1,length(xcm)-1);
    src(nc2:nc2+1) = 1;
    data = [ 1.00 0.99 0.00 0.0  0.01
             1.00 0.01 0.00 0.0  0.99 ];
    mt   = [ ones(1,nc2) ones(1,nc2)*2 ];
    numm = 2; numg = 1; ord = 32; ptype = 0;

else
    disp('NO SUCH PROBLEM')
    return 
end

input   =   struct(       ...
    'numg',         numg, ...     % number of groups FINAL keff = 1.20999183 in 452 iterations
    'numm',         numm, ...     % number of materials
    'xcm',           xcm, ...     % slab bounds
    'xfm',           xfm, ...     % number of fine meshes
    'mt',             mt, ...     % slab material ids
    'data',         data, ...     % mat comp's
    'src',           src, ...     % volume source
    'ord',            ord, ...     % number of ordinates
    'ptype',       ptype, ...     % 0=fixed,1=eigen
    'maxinit',      1000, ...     % max inner iterations
    'maxinerr',     1e-6, ...     % max pointwise phi error
    'maxoutit',     1000, ...     % max iterations    
    'maxouterr',    1e-6, ...     % max pointwise fission density error    
    'maxkerr',      1e-6, ...     % max keff  error1.18739697
    'bcL',             0, ...     % left boundary condition
    'bcR',             0, ...     % right boundary condition
    'acc_e',           0, ...     % eigenvalue acceleration (1=cmfd)FINAL keff = 1.25744374 in 243 iterations
    'acc',             2  ...     % acceleration (0=none,1=cmr,2=cmfd,3=cmsc) 1.1874
);

tic
  [phi,psi,J,x] = sn(input);
toc

hold on
for g = 1:numg
    if g==1
        c = 'g:';
    elseif g==2
        c = 'g:';
    else
        c = 'r:';
    end
    plot(x,phi(:,g),c,'LineWidth',2)
end
plot(x,phi,'k',x,J,'b--','LineWidth',2)
axis([ 0 max(x) min(J)*1.1 max(max(phi))*1.1 ])
grid on

